home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / QFLOAT / QEXP.C < prev    next >
C/C++ Source or Header  |  1996-03-31  |  1KB  |  68 lines
  1. /*                            qexp.c        */
  2. /* exponential function check routine */
  3.  
  4. #include "qhead.h"
  5. #include <stdio.h>
  6. extern QELT qlog2[], qone[];
  7. int qtanh();
  8.  
  9. int qexp( x, y )
  10. QELT *x, *y;
  11. {
  12. QELT num[NQ], den[NQ], x2[NQ];
  13. long i;
  14. QELT sign;
  15.  
  16. /* range reduction theory: x = i + f, 0<=f<1;
  17.  * e**x = e**i * e**f
  18.  * e**i = 2**(i/log 2).
  19.  * Let i/log2 = i1 + f1, 0<=f1<1.
  20.  * Then e**i = 2**i1 * 2**f1, so
  21.  * e**x = 2**i1 * e**(log 2 * f1) * e**f.
  22.  */
  23.  
  24. if( x[1] == 0 )
  25.     {
  26.     qmov( qone, y );
  27.     return(0);
  28.     }
  29. qmov(x, x2);
  30. sign = x2[0];
  31. if( sign != 0 )
  32.     x2[0] = 0;
  33. qifrac( x2, &i, num );        /* x = i + f        */
  34.  
  35. if( i != 0 )
  36.  {
  37.  ltoq( &i, den );        /* floating point i    */
  38.  qdiv( qlog2, den, den );    /* i/log 2        */
  39.  qifrac( den, &i, den );    /* i/log 2  =  i1 + f1    */
  40.  qmul( qlog2, den, den );    /* log 2 * f1        */
  41.  qadd( den, num, x2 );        /* log 2 * f1  + f    */
  42.  }
  43.  
  44. x2[1] -= 1;        /* x/2                */
  45. qtanh( x2, x2 );    /* tanh( x/2 )            */
  46. qadd( x2, qone, num );    /* 1 + tanh            */
  47. qneg( x2 );
  48. qadd( x2, qone, den );    /* 1 - tanh            */
  49. qdiv( den, num, y );    /* (1 + tanh)/(1 - tanh)    */
  50.  
  51. i += y[1];
  52. if( i > MAXEXP )
  53.     {
  54.     if( sign != 0 )
  55.         qclear(y);
  56.     else
  57.         {
  58.         printf( "qexp overflow\n" );
  59.         qinfin(y);
  60.         }
  61.     return 0;
  62.     }
  63. y[1] = i;
  64. if( sign != 0 )
  65.     qdiv( y, qone, y );
  66. return 0;
  67. }
  68.